#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <papi.h>

typedef struct {              // one piece of data in MMFE
  unsigned int size;
  unsigned int* row;
  double* d;
} col;

int P;                    // threads
col* m;                   // matrix
double* v;                // vector
double** r;               // result

unsigned short* iszero;   // index of zeros in vector

// size of the matrix
unsigned int width;
unsigned int height;

void* mult(void* p) {
  unsigned int pid = (unsigned int) p;
  unsigned int i, j, start_col, end_col;
  start_col = width/P * pid;
  end_col = (width/P * (pid+1)) > width ? width : (width/P * (pid+1));

  for(i=start_col; i<end_col; i++) {
    for(j=0; j<m[i].size; j++) {
      r[pid][m[i].row[j]] += m[i].d[j] * v[i];
    }
  }
}

int main(int argc, char* argv[]) {
  FILE *mfile, *vfile;
  unsigned int i, ii, trow, tcol, *t;
  double tdouble;
  char trash[100];
  pthread_t *thread;
  void* status;
  unsigned long start, end;
  unsigned long chkflg = PAPI_VER_CURRENT;

/**************************************************
 *               Init PAPI
**************************************************/
  PAPI_library_init(chkflg);
  if (chkflg != PAPI_VER_CURRENT) {
    printf("Error PAPI Library out of date\n");
    exit(1);
  }

/**************************************************
 *               Vector File
**************************************************/

  // open vfile
  vfile = fopen(argv[2], "r");
  fgets(trash, 100, vfile);
  // get sizes
  fscanf(vfile, "%i %*i %*i", &width);
  // allocate memory
  v = malloc(sizeof(double)*width);
  iszero = malloc(sizeof(unsigned short)*width);
  // init
  for(i=0; i<width; i++) {
    v[i] = 0;
    iszero[i] = 0;
  }

  // populate v;
  while(fscanf(vfile, "%i %*i %lf", &trow, &tdouble)!= EOF) {
    v[trow-1] = tdouble;
  }
  // mark zeros;
  for(i=0; i<width; i++)
    iszero[i] = (v[i] == 0);

  fclose(vfile);

/**************************************************
 *               Matrix File
**************************************************/

  // open mfile
  mfile = fopen(argv[1], "r");
  fgets(trash, 100, vfile);
  // get sizes
  fscanf(mfile, "%i %i %*s", &height, &width);
  // allocate memory
  m = malloc(sizeof(col)*width);
  t = malloc(sizeof(unsigned int)*width);
  // init
  for(i=0; i<width; i++) {
    m[i].size = 0;
    t[i] = 0;
  }

  // count size of each row
  while(fscanf(vfile, "%*i %i %*f", &tcol)!=EOF) {
    if (!iszero[tcol-1]) {
      m[tcol-1].size++;
    }
  }
  // allocate memory
  for(i=0; i<width; i++) {
    m[i].d = malloc(sizeof(double)*m[i].size);
    m[i].row = malloc(sizeof(unsigned int)*m[i].size);
  }

  rewind(vfile);

  // populate m
  fgets(trash, 100, vfile);
  fgets(trash, 100, vfile);
  while(fscanf(vfile, "%i %i %lf", &trow, &tcol, &tdouble)!=EOF) {
    if (!iszero[tcol-1]) {
      m[tcol-1].d[t[tcol-1]] = tdouble;
      m[tcol-1].row[t[tcol-1]] = trow-1;
      t[tcol-1]++;
    }
  }

  free(t);
  fclose(mfile);

//  printf("matrix\n");
//  for(i=0; i<width; i++) {
//    for(j=0; j<m[i].size; j++) {
//      printf("%i  %i  %lf\n", m[i].row[j], i, m[i].d[j]);
//    }
//  }
//  printf("\n");
//
//  printf("vector\n");
//  for(i=0; i<width; i++) {
//    printf("%i %lf\n", i, v[i]);
//  }
//  printf("\n");


/**************************************************
 *              Create Thread
**************************************************/

  // determine thread size
  P = 4;

  // allocate memory
  r = malloc(sizeof(double*)* P);
  for(i=0; i<P; i++) {
    r[i] = malloc(sizeof(double)*height);
    for(ii=0; ii<height; ii++) {
      r[i][ii] = 0;
    }
  }
  thread = malloc(sizeof(pthread_t)*P);

  // create thread
  start = PAPI_get_real_usec();
  for(i=0; i<P; i++) {
    pthread_create(&thread[i], NULL, mult, (void*)i);
  }

  // join threads
  for(i=0; i<P; i++) {
    pthread_join(thread[i], status);
  }

  // combine results
  for(i=0; i<height; i++) {
    //r[0][i] = r[1][i] + r[2][i] + r[3][i];
    r[0][i] += r[1][i];
  }

  end = PAPI_get_real_usec();
  printf("result\n");
  for(i=0; i<height; i++) {
    printf("%lf\n", r[0][i]);
  }
  printf("time: %ld", end-start);

  return 0;
}

